#include "TChain.h"
#include "TTree.h"
#include "TH1F.h"
#include "TH2D.h"
#include "TCanvas.h"
#include "THnSparse.h"

#include "AliAnalysisTask.h"
#include "AliAnalysisManager.h"

#include "AliESDEvent.h"
#include "AliESDInputHandler.h"

#include "AliMCEventHandler.h"
#include "AliMCEvent.h"

#include <AliCentrality.h>

#include "AliESDtrackCuts.h"
#include "AliStack.h"

#include "AliAnalysisTaskLRC.h"

#define PT_MIN 0.3  //set range for Pt
#define PT_MAX 1.5  //(if this selection is performed in ESD track loop)

#define ETA_MIN -0.8  //set range for pseudorapidity
#define ETA_MAX 0.8   //(if this selection is performed in ESD track loop)


ClassImp(AliAnalysisTaskLRC)

//________________________________________________________________________
AliAnalysisTaskLRC::AliAnalysisTaskLRC(const char *name) 
    : AliAnalysisTaskSE(name)
    , fVxMax(0)
    , fVyMax(0)
    , fVzMax(0)
    , fOutput(0)
    , fHistPtPlus(0)
    , fHistPtMinus(0)
    , fHistPt(0)
    , fHistPtPlusAfterPtCut(0)
    , fHistPtMinusAfterPtCut(0)
    , fHistPtAfterPtCut(0)
    , fHistEtaPlus(0)
    , fHistEtaMinus(0)
    , fHistEta(0)
    , fHistEtaPlusAfterPtCut(0)
    , fHistEtaMinusAfterPtCut(0)
    , fHistEtaAfterPtCut(0)
    , fHistEtaAfterCuts(0)
    , fHistPtAfterCuts(0)
    , fHistEtaPtAfterCuts(0)
    , fTrackCuts(0)
    , fHistCharge(0)
    //, fHistY(0)
    //, fHistCent(0)
    , fHistVtxX(0)
    , fHistVtxY(0)
    , fHistVtxZ(0)
    , fHistTPCNcls(0)
    , fHistITSNcls(0)
    , fHistDCAxyz(0)
    , fHistChi2PerClusterTPC(0)
    , fHistNch(0)
    , fHistNchAccepted(0)
    , fHistEtaTPCNcls(0)
    , fHistEtaPhiPlus(0)
    , fHistEtaPhiMinus(0)
    , fHistEtaPhi(0)
    , fHistEtaPhiPlus64(0)
    , fHistEtaPhiMinus64(0)
    , fHistEtaPhi64(0)
    , fHistsSparse(0)
    , fHistNpEvent(0)
    , fHistNmEvent(0)
    //, fHistPtEtaEvent(0)
    , fHistNchEvent(0)
    , fHistNchAcceptedEvent(0)
    , fHistEventCutStats(0)
    , fHistTrackCutStats(0)
{
    // Constructor

    // Define input and output slots here
    // Input slot #0 works with a TChain
    DefineInput(0, TChain::Class());
    // Output slot #0 id reserved by the base class for AOD
    // Output slot #1 writes into a TH1 container
    DefineOutput(1, TList::Class());
}

AliAnalysisTaskLRC::AliAnalysisTaskLRC()
    : AliAnalysisTaskSE("user_task")
    , fVxMax(0)
    , fVyMax(0)
    , fVzMax(0)
    , fOutput(0)
    , fHistPtPlus(0)
    , fHistPtMinus(0)
    , fHistPt(0)
    , fHistPtPlusAfterPtCut(0)
    , fHistPtMinusAfterPtCut(0)
    , fHistPtAfterPtCut(0)
    , fHistEtaPlus(0)
    , fHistEtaMinus(0)
    , fHistEta(0)
    , fHistEtaPlusAfterPtCut(0)
    , fHistEtaMinusAfterPtCut(0)
    , fHistEtaAfterPtCut(0)
    , fHistEtaAfterCuts(0)
    , fHistPtAfterCuts(0)
    , fHistEtaPtAfterCuts(0)
    , fTrackCuts(0)
    , fHistCharge(0)
    //, fHistY(0)
    //, fHistCent(0)
    , fHistVtxX(0)
    , fHistVtxY(0)
    , fHistVtxZ(0)
    , fHistTPCNcls(0)
    , fHistITSNcls(0)
    , fHistDCAxyz(0)
    , fHistChi2PerClusterTPC(0)
    , fHistNch(0)
    , fHistNchAccepted(0)
    , fHistEtaTPCNcls(0)
    , fHistEtaPhiPlus(0)
    , fHistEtaPhiMinus(0)
    , fHistEtaPhi(0)
    , fHistEtaPhiPlus64(0)
    , fHistEtaPhiMinus64(0)
    , fHistEtaPhi64(0)
    , fHistsSparse(0)
    , fHistNpEvent(0)
    , fHistNmEvent(0)
    //, fHistPtEtaEvent(0)
    , fHistNchEvent(0)
    , fHistNchAcceptedEvent(0)
    , fHistEventCutStats(0)
    , fHistTrackCutStats(0)
{
    // Constructor

    // Define input and output slots here
    // Input slot #0 works with a TChain
    DefineInput(0, TChain::Class());
    // Output slot #0 id reserved by the base class for AOD
    // Output slot #1 writes into a TH1 container
    DefineOutput(1, TList::Class());
}

AliAnalysisTaskLRC::~AliAnalysisTaskLRC()
{
    // Destructor
}

//________________________________________________________________________
void AliAnalysisTaskLRC::UserCreateOutputObjects()
{
    // Create histograms
    // Called once
    //
    // For histograms description see .h file

    fOutput = new TList();

    fHistPtPlus = new TH1F("fPtPlus", "+ Particles Vs Pt", 250, 0, 5 );
    fHistPtPlus->GetXaxis()->SetTitle("GeV");
    fHistPtPlus->GetYaxis()->SetTitle("N entries");
    fOutput->Add(fHistPtPlus);

    fHistPtMinus = new TH1F("fPtMinus", "- Particles Vs Pt", 250, 0, 5 );
    fHistPtMinus->GetXaxis()->SetTitle("GeV");
    fHistPtMinus->GetYaxis()->SetTitle("N entries");
    fOutput->Add(fHistPtMinus);

    fHistPt = new TH1F("fPt", "ch Particles Vs Pt", 250, 0, 5 );
    fHistPt->GetXaxis()->SetTitle("GeV");
    fHistPt->GetYaxis()->SetTitle("N entries");
    fOutput->Add(fHistPt);

    fHistPtPlusAfterPtCut = new TH1F("fPtPlusAfterPtCut", "+ Particles Vs Pt after Pt cut", 200, 0, 2 );
    fHistPtPlusAfterPtCut->GetXaxis()->SetTitle("GeV");
    fHistPtPlusAfterPtCut->GetYaxis()->SetTitle("N entries");
    fOutput->Add(fHistPtPlusAfterPtCut);

    fHistPtMinusAfterPtCut = new TH1F("fPtMinusAfterPtCut", "- Particles Vs Pt after Pt cut", 200, 0, 2 );
    fHistPtMinusAfterPtCut->GetXaxis()->SetTitle("GeV");
    fHistPtMinusAfterPtCut->GetYaxis()->SetTitle("N entries");
    fOutput->Add(fHistPtMinusAfterPtCut);

    fHistPtAfterPtCut = new TH1F("fPtAfterPtCut", "ch Particles Vs Pt after Pt cut", 200, 0, 2 );
    fHistPtAfterPtCut->GetXaxis()->SetTitle("GeV");
    fHistPtAfterPtCut->GetYaxis()->SetTitle("N entries");
    fOutput->Add(fHistPtAfterPtCut);

    fHistEtaPlus = new TH1F("fEtaPlus", "+ particles vs Eta", 300, -1.5, 1.5 );
    fHistEtaPlus->GetXaxis()->SetTitle("#eta");
    fHistEtaPlus->GetYaxis()->SetTitle("N entries");
    fOutput->Add(fHistEtaPlus);

    fHistEtaMinus = new TH1F("fEtaMinus", "- particles vs Eta", 300, -1.5, 1.5 );
    fHistEtaMinus->GetXaxis()->SetTitle("#eta");
    fHistEtaMinus->GetYaxis()->SetTitle("N entries");
    fOutput->Add(fHistEtaMinus);

    fHistEta = new TH1F("fEta", "ch particles vs Eta", 300, -1.5, 1.5 );
    fHistEta->GetXaxis()->SetTitle("#eta");
    fHistEta->GetYaxis()->SetTitle("N entries");
    fOutput->Add(fHistEta);

    fHistEtaPlusAfterPtCut = new TH1F("fEtaPlusAfterPtCut", "+ particles vs Eta after Pt Cut", 300, -1.5, 1.5 );
    fHistEtaPlusAfterPtCut->GetXaxis()->SetTitle("#eta");
    fHistEtaPlusAfterPtCut->GetYaxis()->SetTitle("N entries");
    fOutput->Add(fHistEtaPlusAfterPtCut);

    fHistEtaMinusAfterPtCut = new TH1F("fEtaMinusAfterPtCut", "- particles vs Eta after Pt Cut", 300, -1.5, 1.5 );
    fHistEtaMinusAfterPtCut->GetXaxis()->SetTitle("#eta");
    fHistEtaMinusAfterPtCut->GetYaxis()->SetTitle("N entries");
    fOutput->Add(fHistEtaMinusAfterPtCut);

    fHistEtaAfterPtCut = new TH1F("fEtaAfterPtCut", "ch particles vs Eta after Pt Cut", 300, -1.5, 1.5 );
    fHistEtaAfterPtCut->GetXaxis()->SetTitle("#eta");
    fHistEtaAfterPtCut->GetYaxis()->SetTitle("N entries");
    fOutput->Add(fHistEtaAfterPtCut);

    fHistEtaAfterCuts = new TH1F("fEtaAfterCuts", "ch particles vs Eta after all cuts", 300, -1.5, 1.5 );
    fHistEtaAfterCuts->GetXaxis()->SetTitle("#eta");
    fHistEtaAfterCuts->GetYaxis()->SetTitle("N entries");
    fOutput->Add(fHistEtaAfterCuts);

    fHistPtAfterCuts = new TH1F("fPtAfterCuts", "ch Particles Vs Pt after all cuts", 200, 0, 2 );
    fHistPtAfterCuts->GetXaxis()->SetTitle("GeV");
    fHistPtAfterCuts->GetYaxis()->SetTitle("N entries");
    fOutput->Add(fHistPtAfterCuts);

    fHistEtaPtAfterCuts = new TH2F("fEtaPtAfterCuts", "ch Particles Vs Eta and Pt after all cuts", 200, -1, 1, 200, 0, 2);
    fHistEtaPtAfterCuts->GetXaxis()->SetTitle("#eta");
    fHistEtaPtAfterCuts->GetYaxis()->SetTitle("GeV");
    fOutput->Add(fHistEtaPtAfterCuts);

    fHistCharge = new TH1I("fCharge", "N particles vs charge", 5, -2.5, 2.5 );
    fHistCharge->GetXaxis()->SetTitle("charge");
    fHistCharge->GetYaxis()->SetTitle("N entries");
    fOutput->Add(fHistCharge);

    fHistVtxX = new TH1F("fVtxX", "Vertex x coordinate distribution", 100, -0.5, 0.5 );
    fHistVtxX->GetXaxis()->SetTitle("Vx (cm)");
    fHistVtxX->GetYaxis()->SetTitle("N entries");
    fOutput->Add(fHistVtxX);

    fHistVtxY = new TH1F("fVtxY", "Vertex y coordinate distribution", 100, -0.5, 0.5 );
    fHistVtxY->GetXaxis()->SetTitle("Vy (cm)");
    fHistVtxY->GetYaxis()->SetTitle("N entries");
    fOutput->Add(fHistVtxY);

    fHistVtxZ = new TH1F("fVtxZ", "Vertex z coordinate distribution", 200, -10, 10 );
    fHistVtxZ->GetXaxis()->SetTitle("Vz (cm)");
    fHistVtxZ->GetYaxis()->SetTitle("N entries");
    fOutput->Add(fHistVtxZ);

    fHistTPCNcls = new TH1I("fTPCNcls", "N of TPC clusters distribution", 200, 0, 200 );
    fHistTPCNcls->GetXaxis()->SetTitle("N of TPC clusters after cuts");
    fHistTPCNcls->GetYaxis()->SetTitle("N entries");
    fOutput->Add(fHistTPCNcls);

    fHistITSNcls = new TH1I("fITSNcls", "N of ITS clusters distribution", 12, -2, 10 );
    fHistITSNcls->GetXaxis()->SetTitle("N of ITS clusters after cuts");
    fHistITSNcls->GetYaxis()->SetTitle("N entries");
    fOutput->Add(fHistITSNcls);

    fHistDCAxyz = new TH2F("fDCAxyz", "DCAxy vs DCAz", 28, -0.7, 0.7, 28, -0.7, 0.7);
    fHistDCAxyz->GetXaxis()->SetTitle("DCA to vertex XY");
    fHistDCAxyz->GetYaxis()->SetTitle("DCA to Vertex Z");
    fOutput->Add(fHistDCAxyz);

    fHistChi2PerClusterTPC = new TH1F("fChi2PerClusterTPC", "Chi^2 per cluster TPC", 60, 0, 6);
    fHistChi2PerClusterTPC->GetXaxis()->SetTitle("Chi^2/Nclusters TPC");
    fOutput->Add(fHistChi2PerClusterTPC);

    fHistNch = new TH1F("fNch","N of Nch particles in event", 200, 0, 200);
    fHistNch->GetXaxis()->SetTitle("N of charged tracks");
    fHistNch->GetYaxis()->SetTitle("N of events");
    fOutput->Add(fHistNch);

    fHistNchAccepted = new TH1F("fNchAccepted","N of accepted Nch particles in event", 100, 0, 100);
    fHistNchAccepted->GetXaxis()->SetTitle("N of charged tracks accepted");
    fHistNchAccepted->GetYaxis()->SetTitle("N of events");
    fOutput->Add(fHistNchAccepted);

    fHistEtaTPCNcls = new TH2F("fEtaTPCNcls","N of charged particles Vs Eta and TPCNcls", 60, -1.5, 1.5, 100, 70, 170);
    fHistEtaTPCNcls->GetXaxis()->SetTitle("#eta");
    fHistEtaTPCNcls->GetYaxis()->SetTitle("N of TPC clusters");
    fOutput->Add(fHistEtaTPCNcls);

    fHistEtaPhiPlus = new TH2F("fEtaPhiPlus", "+ particles vs EtaPhi after all cuts", 20, -1, 1, 32, 0, 2*TMath::Pi() );
    fHistEtaPhiPlus->GetXaxis()->SetTitle("#eta");
    fHistEtaPhiPlus->GetYaxis()->SetTitle("#phi");
    fOutput->Add(fHistEtaPhiPlus);

    fHistEtaPhiMinus = new TH2F("fEtaPhiMinus", "- particles vs EtaPhi after all cuts", 20, -1, 1, 32, 0, 2*TMath::Pi() );
    fHistEtaPhiMinus->GetXaxis()->SetTitle("#eta");
    fHistEtaPhiMinus->GetYaxis()->SetTitle("#phi");
    fOutput->Add(fHistEtaPhiMinus);

    fHistEtaPhi = new TH2F("fEtaPhi", "N particles vs EtaPhi after all cuts", 20, -1, 1, 32, 0, 2*TMath::Pi() );
    fHistEtaPhi->GetXaxis()->SetTitle("#eta");
    fHistEtaPhi->GetYaxis()->SetTitle("#phi");
    fOutput->Add(fHistEtaPhi);

    fHistEtaPhiPlus64 = new TH2F("fEtaPhiPlus64", "+ particles vs EtaPhi after all cuts", 10, -1, 1, 8, 0, 2*TMath::Pi() );
    fHistEtaPhiPlus64->GetXaxis()->SetTitle("#eta");
    fHistEtaPhiPlus64->GetYaxis()->SetTitle("#phi");
    fOutput->Add(fHistEtaPhiPlus64);

    fHistEtaPhiMinus64 = new TH2F("fEtaPhiMinus64", "- particles vs EtaPhi after all cuts", 10, -1, 1, 8, 0, 2*TMath::Pi() );
    fHistEtaPhiMinus64->GetXaxis()->SetTitle("#eta");
    fHistEtaPhiMinus64->GetYaxis()->SetTitle("#phi");
    fOutput->Add(fHistEtaPhiMinus64);

    fHistEtaPhi64 = new TH2F("fEtaPhi64", "N particles vs EtaPhi after all cuts", 10, -1, 1, 8, 0, 2*TMath::Pi() );
    fHistEtaPhi64->GetXaxis()->SetTitle("#eta");
    fHistEtaPhi64->GetYaxis()->SetTitle("#phi");
    fOutput->Add(fHistEtaPhi64);

    /*fHistY = new TH1F("fY", "N particles vs Y", 400, -2, 2 );
  fHistY->GetXaxis()->SetTitle("Y");
  fHistY->GetYaxis()->SetTitle("N entries");
  fOutput->Add(fHistY);*/

    /*fHistCent = new TH1F("fCent", "N events vs centrality class", 110, -10, 100);
  fHistCent->GetXaxis()->SetTitle("centrality class");
  fHistCent->GetYaxis()->SetTitle("N events");
  fOutput->Add(fHistCent);*/

    //fOutput->Add(fTrackCuts);

    //for PbPb:
    /*Int_t bins[11]={100,800,800,800,800,1600,1600,400,400,140,140};//number of bins
  Double_t min[11]={0,-0.5,-0.5,-0.5,-0.5,-0.5,-0.5,-200.5,-200.5,0.195,0.195};//min values
  Double_t max[11]={100,799.5,799.5,799.5,799.5,1599.5,1599.5,199.5,199.5,1.595,1.595};//max values*/

    //for pp:       {Np_B,Np_F,Nm_B,Nm_F,Nch_B,Nch_F,Q_B, Q_F,Q_F*Q_B,Q_F*Q_F }

    TString AxisNames[10] = {"Np_B","Np_F","Nm_B","Nm_F","Nch_B","Nch_F","Q_B", "Q_F","Q_F*Q_B","Q_F*Q_F"};
    Int_t bins[10] = {  30,  30,  30,  30, 60, 60,   50,   50,   300,   300};//number of bins
    Double_t min[10]={-0.5,-0.5,-0.5,-0.5,-0.5,-0.5,-25.5,-25.5, -150.5, -0.5};//min values
    Double_t max[10]={29.5,29.5,29.5,29.5,59.5,59.5, 24.5, 24.5,  149.5, 299.5};//max values


    fHistsSparse = new THnSparseF*[20];


        for (Int_t i=0;i<=15;i++)
        {

            TString strSparse = Form( "fPhi2piEta%d", i+1);
            fHistsSparse[i] = new THnSparseF(strSparse,"Sparse hist",10,bins,min,max);

            for (Int_t d=0; d<=9;d++)
            {
                fHistsSparse[i]->GetAxis(d)->SetNameTitle(AxisNames[d],AxisNames[d]);

            }

            fOutput->Add(fHistsSparse[i]);

        }


    fHistNpEvent = new TH1F("fNpEvent", "Nplus in one event vs EtaPhi", 18, -0.9, 0.9);

    fHistNmEvent = new TH1F("fNmEvent", "Nminus in one event vs EtaPhi", 18, -0.9, 0.9);

    //fHistPtEtaEvent = new TH2F("f2PtEtaEvent", "Pt vs Eta in one event", 140, 0.195, 1.595, 16, -0.8, 0.8);

    fHistNchEvent = new TH1I("fNchEvent", "Counter of Ncharged in one event ", 1, 0.5, 1.5 );

    fHistNchAcceptedEvent = new TH1I("fNchAcceptedEvent", "Counter of Ncharged accepted in one event ", 1, 0.5, 1.5 );

    //Event level QA
    Int_t nStatBins = 8;
    fHistEventCutStats = new TH1D("fHistEventCutStats","Event statistics;;N_{events}", nStatBins,-0.5,nStatBins-0.5);
    TString gEventCutBinNames[8] = {"Total","No trigger","Centrality","No vertex","Bad vertex position","fESD not available","LowMult cut","Analyzed"};

    for(Int_t i = 1; i <= nStatBins; i++)
    {
        fHistEventCutStats->GetXaxis()->SetBinLabel(i,gEventCutBinNames[i-1].Data());
    }
    fOutput->Add(fHistEventCutStats);

    //Track level QA
    fHistTrackCutStats = new TH1D("fHistTrackCutStats","Event statistics;;N_{events}", 7,-0.5,6.5);
    TString gTrackCutBinNames[7] = {"Total","AliESDtrackCuts","HighPt cut","LowPt cut","High Eta cut","Low Eta Cut","Analyzed"};

    for(Int_t i = 1; i <= 7; i++)
    {
        fHistTrackCutStats->GetXaxis()->SetBinLabel(i,gTrackCutBinNames[i-1].Data());
    }
    fOutput->Add(fHistTrackCutStats);

    //fOutput->Add(fTrackCuts);
}

//________________________________________________________________________
void AliAnalysisTaskLRC::UserExec(Option_t *) 
{
    // Main loop
    // Called for each event

    //if ( fHaveMC )
    //	UserExecMC();
    AliESDEvent *fESD = dynamic_cast<AliESDEvent*>(InputEvent());
    fHistEventCutStats->Fill(0); //total number of events

    if ( !fESD )
    {
        Printf("ERROR: fESD not available\n");
        fHistEventCutStats->Fill(5);
        return;
    }


    //Trigger
    Bool_t lTrigger = kTRUE;
    if( fESD )
    {
        lTrigger = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
        if( !lTrigger )
        {
            fHistEventCutStats->Fill(1);
            PostData(1, fOutput);
            return;
        }
    }

    // Vertex present

    const AliESDVertex *vertex = fESD->GetPrimaryVertex();

    Bool_t lVertexPresent = ( (vertex) && (vertex->GetNContributors() > 0 ) && ( vertex->GetZRes() != 0 ) );
    if( !lVertexPresent )
    {
        fHistEventCutStats->Fill(3); //no vertex
        PostData(1, fOutput);
        return;
    }


    // Vertex in range
    Bool_t lVertexAcceptable = (TMath::Abs(vertex->GetXv()) < fVxMax) && (TMath::Abs(vertex->GetYv()) < fVyMax);
    if( lVertexAcceptable )
    {
        if( fVzMax>0 )
            lVertexAcceptable = (TMath::Abs(vertex->GetZv()) < fVzMax);
        else
            lVertexAcceptable = (TMath::Abs(vertex->GetZv()) > -fVzMax);
    }

    if(!lVertexAcceptable)
    {
        fHistEventCutStats->Fill(4); //bad vertex
        PostData(1, fOutput);
        return;
    }

    Double_t lVtxX = vertex->GetXv();
    fHistVtxX->Fill(lVtxX);//x coordinate distribution of vertex

    Double_t lVtxY = vertex->GetYv();
    fHistVtxY->Fill(lVtxY);//y coordinate distribution of vertex

    Double_t lVtxZ = vertex->GetZv();
    fHistVtxZ->Fill(lVtxZ);//z coordinate distribution of vertex


    //Printf("There are %d tracks in this event\n", fESD->GetNumberOfTracks());

    //centrality selection (only for AA)
    /*AliCentrality *centrality = fESD->GetCentrality();
   Double_t lCentrality = centrality->GetCentralityPercentile("V0M");

   // GetCentralityClass10 returns the number of 10% centrality class, an integer from 0 to 9
   // "V0M" means that we using V0 detector as the default estimator for centrality
   // other possible methods:
   // GetCentralityClass5 - number of 5% centrality class
   // GetCentralityPercentile - percentile, a float from 0 to 100 (or to the trigger efficiency)

  if (lCentrality < 0.25 || lCentrality > 80) //select centrality if needed
  {
  fHistEventCutStats->Fill(2);
  return;
  }

  //centrality classes if needs:
  Int_t lCentralityClass; //centrality class from 0-5,5-10,10-20 to 70-80

  if (lCentralityClass10 > 0 && lCentralityClass10 < 8)
  {lCentralityClass = lCentralityClass10 + 1;}
  else if (lCentralityClass10 == 0)
  {lCentralityClass = centrality->GetCentralityClass5("V0M");}*/


    // user trigger to be added (if needed)
    // user vertex cut (if needed)

    fHistNpEvent->Reset();
    fHistNmEvent->Reset();
    //fHistPtEtaEvent->Reset();
    fHistNchEvent->Reset();
    fHistNchAcceptedEvent->Reset();


    // ESD Track loop
    for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++)
    {
        AliESDtrack* track = fESD->GetTrack(iTracks);

        fHistTrackCutStats->Fill(0); //total number of raw tracks

        if( !fTrackCuts->AcceptTrack(track) )
        {
            fHistTrackCutStats->Fill(1); //rejectred by cuts
            continue;
        }

        if ( !track )
        {
            printf("ERROR: Could not receive track %d\n", iTracks);
            continue;
        }


        Int_t NclustersTPC = track->GetTPCNcls();//number of TPC clusters distribution AFTER cuts
        fHistTPCNcls->Fill(NclustersTPC);

        Int_t NclustersITS = track->GetITSclusters(0);//number of ITS clusters distribution AFTER cuts
        fHistITSNcls->Fill(NclustersITS);


        // getting the track to vertex parameters
        Float_t b[2];
        Float_t bCov[3];
        track->GetImpactParameters(b,bCov);
        Float_t dcaToVertexXY = b[0];
        Float_t dcaToVertexZ = b[1];
        fHistDCAxyz->Fill(dcaToVertexXY,dcaToVertexZ);

        //chi2 per cluster TPC
        if (NclustersTPC!=0)
        {
            Float_t  chi2PerClusterTPC = track->GetTPCchi2()/(NclustersTPC);
            fHistChi2PerClusterTPC->Fill(chi2PerClusterTPC);
        }

        fHistNchEvent->Fill(1); //count total number of tracks

        Double_t lTrackPt = track->Pt();
        Double_t lTrackEta = track->Eta();
        Int_t lTrackCharge = track->Charge();


        if (lTrackCharge>0)
        {
            fHistPtPlus->Fill(lTrackPt);
            fHistEtaPlus->Fill(lTrackEta);
        }

        if (lTrackCharge<0)
        {
            fHistPtMinus->Fill(lTrackPt);
            fHistEtaMinus->Fill(lTrackEta);
        }

        fHistPt->Fill(lTrackPt);
        fHistEta->Fill(lTrackEta);
        fHistEtaTPCNcls->Fill(lTrackEta, NclustersTPC);

        if ( lTrackPt > PT_MAX)
        {fHistTrackCutStats->Fill(2);}
        if ( lTrackPt < PT_MIN)
        {fHistTrackCutStats->Fill(3);}

        //cut on Pt
        if( lTrackPt < PT_MIN || lTrackPt > PT_MAX )
            continue;

        fHistPtAfterPtCut->Fill(lTrackPt);
        fHistEtaAfterPtCut->Fill(lTrackEta);

        if (lTrackCharge > 0)
        {
            fHistPtPlusAfterPtCut->Fill(lTrackPt);
            fHistEtaPlusAfterPtCut->Fill(lTrackEta);
        }

        if (lTrackCharge < 0)
        {
            fHistPtMinusAfterPtCut->Fill(lTrackPt);
            fHistEtaMinusAfterPtCut->Fill(lTrackEta);
        }


        if ( lTrackEta > ETA_MAX)
        {fHistTrackCutStats->Fill(4);}
        if ( lTrackEta < ETA_MIN)
        {fHistTrackCutStats->Fill(5);}

        //cut on Eta
        if( lTrackEta < ETA_MIN || lTrackEta > ETA_MAX )
            continue;


        fHistEtaAfterCuts->Fill(lTrackEta);
        fHistPtAfterCuts->Fill(lTrackPt);
        fHistEtaPtAfterCuts->Fill(lTrackEta, lTrackPt);

        //Double_t lTrackY = track->Y();
        Double_t lTrackPhi = track->Phi();

        if (lTrackCharge > 0) {fHistNpEvent->Fill(lTrackEta);}
        if (lTrackCharge < 0) {fHistNmEvent->Fill(lTrackEta);}

        //fHistPtEtaEvent->Fill(lTrackPt,lTrackEta);

        fHistNchAcceptedEvent->Fill(1); //count total number of accepted tracks

        fHistEtaPhi->Fill( lTrackEta, lTrackPhi);
        if (lTrackCharge > 0) { fHistEtaPhiPlus->Fill( lTrackEta, lTrackPhi);}
        if (lTrackCharge < 0) { fHistEtaPhiMinus->Fill( lTrackEta, lTrackPhi);}

        fHistEtaPhi64->Fill( lTrackEta, lTrackPhi);
        if (lTrackCharge > 0) { fHistEtaPhiPlus64->Fill( lTrackEta, lTrackPhi);}
        if (lTrackCharge < 0) { fHistEtaPhiMinus64->Fill( lTrackEta, lTrackPhi);}

        fHistCharge->Fill(lTrackCharge);
        //fHistY->Fill( lTrackY );
        fHistTrackCutStats->Fill(6);
    }

    Double_t NchAcceptedEvent = fHistNchAcceptedEvent->Integral(1,1);
    fHistNchAccepted->Fill(NchAcceptedEvent);

    Double_t NchEvent = fHistNchEvent->Integral(1,1);
    fHistNch->Fill(NchEvent);

    //arrays of bin numbers for 16 eta windows
    Int_t downB[16]={ 8, 7, 6, 5, 4, 3, 2,  6, 5, 4, 3, 2,  4, 3, 2,  2};
    Int_t upB[16] = { 9, 8, 7, 6, 5, 4, 3,  9, 8, 7, 6, 5,  9, 8, 7,  9};
    Int_t downF[16]={ 10,11,12,13,14,15,16, 10,11,12,13,14, 10,11,12, 10};
    Int_t upF[16] = { 11,12,13,14,15,16,17, 13,14,15,16,17, 15,16,17, 17};


        for (Int_t i=0;i<=15;i++)
        {
            Double_t Np_B = fHistNpEvent->Integral(downB[i],upB[i]);//Np backward
            Double_t Np_F = fHistNpEvent->Integral(downF[i],upF[i]);//Np forward
            Double_t Nm_B = fHistNmEvent->Integral(downB[i],upB[i]);//Nm backward
            Double_t Nm_F = fHistNmEvent->Integral(downF[i],upF[i]);//Nm forward
            Double_t Nch_B = Np_B + Nm_B;
            Double_t Nch_F = Np_F + Nm_F;
            Double_t Q_B = Np_B - Nm_B;
            Double_t Q_F = Np_F - Nm_F;
            //Double_t PtMean_B = fHistPtEtaEvent->ProjectionX("name2",downB[i],upB[i])->GetMean();//backward
            //Double_t PtMean_F = fHistPtEtaEvent->ProjectionX("name",downF[i],upF[i])->GetMean();//forward
            Double_t x[10]={Np_B, Np_F, Nm_B, Nm_F, Nch_B, Nch_F, Q_B, Q_F, Q_B*Q_F, Q_F*Q_F};
            fHistsSparse[i]->Fill(x);
        }


    //fHistCent->Fill(lCentrality);
    fHistEventCutStats->Fill(7); //analyzed

    PostData(1, fOutput);
}



//________________________________________________________________________
void AliAnalysisTaskLRC::Terminate(Option_t *) 
{
    // Draw result to the screen
    // Called once at the end of the query

    fOutput = dynamic_cast<TList*> (GetOutputData(1));
    if (!fOutput) {
        //printf("ERROR: Output list not available\n");
        return;
    }

}
